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ABSTRACT 

In  this  paper  we  study  uniform  high  order  spectral  methods  to  solve  multi-dimensional 
Euler  equations  for  gas  dynamics.  Uniform  high  order  spectral  approximations  with  spectral 
accuracy  in  smooth  regions  of  solutions  are  constructed  by  introducing  the  idea  of  the  Essen¬ 
tially  Non-Oscillatory  (ENO)  polynomial  interpolations  into  the  spectral  methods.  Based 
on  the  new  approximations,  we  propose  nonoscillatory  spectral  methods  which  possess  the 
properties  of  both  upwinding  difference  schemes  and  spectral  methods.  We  present  numerical 
results  for  the  inviscid  Burgers’  equation,  and  for  one  dimensional  Euler  equations  including 
the  interactions  between  a  shock  wave  and  density  disturbance,  Sod’s  and  Lax’s  shock  tube 
problems,  and  the  blast  wave  problem.  Finally,  we  simulate  the  interaction  between  a  Mach 
3  two  dimensional  shock  wave  and  a  rotating  vortex. 
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1  Introduction 


Recently,  high  order  numerical  methods  have  attracted  considerable  attention^  for  the  sim- 

) 

ulation  of  flows  with  shock  waves  and  different  scales,  especially  for  turbulent  flows  affected 
by  shock  wave  interactions.  Those  high  order  methods  are  expected  to  produce  nonoscil- 
latory  sharp  shock  profiles  without  excessive  overall  numerical  diffusion  and,  at  the  same 
time,  be  able  to  resolve  the  small  scales  of  the  flow  field  elsewhere.  Recent  results  with 
essentially  nonoscillatory  (ENO)  finite  difference  methods  have  shown  considerable  progress 
in  this  direction  [9],  [17].  Spectral  methods,  as  high  order  global  methods,  have  been  very 
successful  in  the  study  of  turbulent  flows  and  flow  transition  problems  when  the  solutions  of 
the  fluid  problems  are  smooth.  For  those  problems,  spectral  methods  have  been  shown  to 
have  an  accuracy  higher  than  any  algebraic  order  (so  called  spectral  accuracy)  [5].  However 
it  remains  to  show  that  spectral  methods  will  also  be  successful  in  computing  flows  with 
shock  waves. 

In  this  paper,  we  continue  our  previous  work  [4]  in  designing  essentially  nonoscillatory 
spectral  methods  for  computing  the  weak  solutions  of  the  hyperbolic  system  of  conservation 
laws 

ut  +  f(u),  +  g(u)w  =  0  (1-1) 

u(x,y,0)  =  Uo(z,s/).  (1.2) 

Here,  as  usual,  u  =  (uj,  •  •  •  ,ua)T  is  a  state  vector  and  f(u),g(u)  are  the  vector- valued 
flux  functions  of  s  components.  The  system  is  assumed  to  be  hyperbolic  in  the  sense  that 
for  any  real  vector  (  =  (^x ,  ^2),  the  matrix  £ifj;  +  always  has  s  -  real  eigenvalues  and 
a  complete  set  of  eigenvectors.  The  solutions  to  (1.1)  usually  develop  discontinuities  in  the 
form  of  shock  waves  and  contact  discontinuities. 

In  applying  spectral  methods  to  problems  having  discontinuous  solutions,  a  key  issue  is 
how  to  deal  with  the  Gibbs  phenomenon  caused  by  the  discontinuities  of  the  solutions.  The 
overall  accuracy  of  spectral  methods  will  be,  at  most,  first  order  everywhere  in  the  presence 
of  Gibbs  oscillations.  There  are  various  filtering  techniques  to  recover  spectral  accuracy  in 
regions  away  from  the  discontinuities  [8],  [14].  On  the  other  hand,  one-sided  filtering  can 
be  used  to  obtain  uniform  convergence  in  regions  close  to  the  discontinuities^].  As  another 
approach  to  treat  the  Gibbs  oscillation,  in  [4]  we  proposed  a  nonoscillatory  spectral  approxi¬ 
mation  to  discontinuous  solutions  by  adding  piecewise  linear  functions,  such  as  sawtooth-like 
functions  and  step  functions,  to  the  conventional  Fourier  trigonometric  or  Chebyshev  poly¬ 
nomial  spaces.  Those  additional  functions  are  usuu  to  resolve  the  discontinuities  in  the 


solutions  caused  by  shock  waves  and  contact  discontinuities.  The  cell-averaged  form  of 
(1.1)  is  used  to  formulate  the  numerical  schemes,  resulting  in  Godunov-type  shock  capturing 
algorithms.  The  usual  reconstruction  step  between  cell  averages  and  point  values  of  the  nu¬ 
merical  solutions  in  such  schemes  can  be  performed  efficiently  with  Fast  Fourier  Transforms. 
However,  a  common  problem  with  cell-averaged  formulation  is  the  costly  implementation  of 
the  reconstruction  in  multi-dimensional  problems. 

In  this  paper,  we  adopt  the  same  philosophy  as  in  [4],  however,  a  more  robust  and 
sophisticated  technique  will  be  introduced.  With  the  new  technique,  we  will  be  able  to 
achieve  global  convergence  up  to  any  given  m-th  order  (m  >  0)  and  ,  meanwhile,  retain 
spectral  accuracy  in  the  regions  away  from  the  discontinuities.  In  order  to  achieve  these  goals, 
we  incorporate  the  main  idea  of  the  ENO  polynomial  interpolations  [9]  into  our  construction 
of  uniform  spectral  approximations.  We  also  introduce  the  idea  of  upwind  differencing  from 
conservative  finite  difference  methods  into  the  design  of  the  spectral  schemes.  The  idea 
of  upwinding  has  proven  very  successful  in  capturing  shock  wave  fronts  and  producing  the 
right  entropy  solutions.  By  using  local  Riemann  solvers  and  flux  limiters,  modern  shock 
capturing  finite  difference  schemes,  such  as  the  TVD  schemes  [9]  ,  MUSCL  type  schemes 
[19],  FCT  schemes  [1],  and  the  more  rerent  ENO  schemes  [9]  [17],  produce  very  satisfactory 
shock  profiles  and  entropy  satisfying  solutions.  The  nonoscillatory  spectral  approximations 
proposed  in  this  paper  will  enable  us  to  bring  the  upwinding  idea  into  the  framework  of 
spectral  methods.  Meanwhile,  the  spectral  schemes  will  be  based  directly  on  the  conservation 
laws  (1-1),  not  on  its  cell-averaged  form.  Thus,  generalization  to  multi-dimensional  cases  is 
straightforward. 

For  systems  of  conservation  laws,  in  order  to  achieve  sharp  shock  profile  without  spurious 
oscillations,  numerical  flux  operators  for  the  scalar  equations  are  usually  applied  to  the 
locally  defined  characteristic  variables.  Because  of  this  complication,  it  has  been  realized 
that  the  cost  of  upwinding  schemes  is  much  greater  than  that  of  the  centered  difference 
schemes.  Several  attempts  have  been  made  to  eliminate  this  shortcoming  by  combining 
centered  difference  schemes  and  upwinding  schemes.  In  [13],  a  mixed  method  of  centered 
difference  schemes  and  ENO  schemes  was  studied  and,  in  [6],  the  authors  suggested  a  type 
of  nonlinear  filtering  technique  to  modify  the  results  of  the  Lax-Wendroff  scheme  at  each 
time  step  to  produce  nonoscillatory  TVD  solutions.  Even  though  it  is  achieved  in  a  quite 
different  way  from  those  in  [6]  and  [13],  the  result  in  this  paper  will  provide  another  example 
of  blending  the  nice  properties  of  both  upwinding  schemes  and  centered  difference  schemes 
(in  this  case,  spectral  schemes). 

This  paper  is  organized  as  follows:  in  Section  2,  we  first  briefly  review  the  method 
proposed  in  [4],  then  present  the  new  method  of  constructing  uniform  convergent,  up  to  any 
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given  m-th  order(7n  >  0),  spectral  approximations  to  discontinuous  functions.  In  Section 

3,  we  study  the  nonoscillatory  spectral  methods  for  scalar  conservation  laws.  Extensions  to 
the  system  of  conservation  laws  and  multi- dimensional  problems  will  be  discussed  in  Section 

4.  In  Section  5,  we  present  numerical  experiments  for  the  new  methods.  First,  the  uniform 
convergence  and  the  spectral  accuracy  of  the  proposed  spectral  approximations  are  tested 
on  discontinuous  functions.  Then  we  study  the  global  accuracy  of  the  spectral  schemes  on 
a  scalar  inviscid  Burgers’  equation  and  on  one  dimensional  Euler  equations  which  model 
the  interaction  of  a  pure  shock  wave  with  density  waves.  Also  we  apply  the  scheme  to 
the  standard  Sod’s  and  Lax’s  test  problems  [16]  in  order  to  check  the  convergence  of  our 
spectral  schemes  to  the  correct  entropy  solutions.  High  order  numerical  results  will  also  be 
presented  for  solving  the  interaction  between  two  blast  waves  [20].  Finally  we  apply  the 
spectral  schemes  to  simulate  the  interactions  between  a  Mach  3  two-dimensional  shock  wave 
and  a  rotating  vortex. 

2  Uniform  High  Order  Spectral  Approximations 

The  conventional  Fourier  spectral  space  has  basis  functions  {e,fcx}|fc|<w.  The  Fourier  expan¬ 
sions  for  discontinuous  functions  converge  very  slowly.  For  instance,  consider  a  sawtooth-like 
function 


F{x,*.,A)  =  a\-, 

V 


for  x  <  x,, 
for  x  >  x,, 


(2.1) 


where  x,  is  the  location  of  the  discontinuity  and  A  —  F^x‘  ^  =  [^]x,  is  the  jump  of 
F(x,x,,A)  across  x,.  The  partial  sum  of  the  Fourier  expansion  of  F(x,xt,A)  is 


Fn(x,x„A)=  £  fk(x.,A)eik*,  (2.2) 

\k\<N 

where 

=  sf  A)e~'k‘dX  =  A{  (Xx.)  for  jfL-o.1,  <2-3> 

From  (2.3)  we  can  see  that  the  Fourier  coefficients  fk(x,,A)  only  decay  like  0(£)  as 
k  — *  oo.  As  a  result,  the  convergence  of  (2.2)  will  be  only  first  order,  and  moreover,  the 
Gibbs  oscillations  near  x,  will  be  in  the  order  of  0(1).  In  order  to  get  rid  of  the  Gibbs 
oscillations,  in  [4]  we  proposed  a  technique  to  construct  essentially  nonoscillatory  spectral 
approximations,  which  we  review  below. 

Let  tt(x)  be  a  piecewise  C°°  periodic  function  with  a  jump  discontinuity  at  x,  with  jump 
[u]x,  and  let  un(x)  be  its  finite  Fourier  expansion,  the  nonoscillatory  spectral  approximation 
is  defined  by 
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(2.4) 


“wM  =  E  +  E  4e"“'e'*1' 

|fc|<W  |fe|>W  lK 

where  y  is  an  approximation  of  x,  and  A1  is  an  approximation  of  [ulT,  and 

f  u(x)e~'kx  dx. 

2k  Jo 

Since  the  second  sum  in  (2.4)  is  actually  F(x,y,A')  —  FN(x,y,  A'),  we  have 

««,(*)=  E  K  -  A(»,  /l')]eifc*  +  F(x,y,A‘).  (2.5) 

|Jb|<JV 

Therefore  u*N(x)  defines  an  approximation  in  the  spectral  space  {e‘fcx}|*.|<jv  augmented 
by  sawtooth-like  functions  F(x,y,A'). 

The  approximation  defined  in  (2.5)  yields  nonoscillatory  numerical  results  for  discon¬ 
tinuous  functions,  and  spectral  schemes  using  this  approximation  have  given  high  order 
accuracy  for  one-dimensional  Euler  gas  dynamics  equations  ([2],  [3]).  In  order  for  (2.5)  to 
be  nonoscillatory,  the  approximation  for  the  location  and  the  magnitude  of  the  shock  should 
be  reasonably  accurate.  Second  order  accuracy  in  the  location  and  first  order  accuracy  in 
the  magnitude  are  needed  to  ensure  the  uniform  nonoscillatory  convergence. 

In  what  follows,  we  present  a  different  method  which  will  be  uniformly  convergent  up  to 
any  given  order  m  >  0  and  ,  at  the  same  time,  retain  the  spectral  accuracy  in  the  smooth 
regions  away  from  the  discontinuities.  Furthermore,  the  requirement  of  accuracy  in  shock 
locations  will  be  much  relaxed  and  an  approximation  to  the  magnitude  of  the  shock  is  not 
needed.  This  makes  the  computation  more  robust. 

Before  we  discuss  the  new  approximation  method  we  review  two  techniques  to  be  used 
in  our  construction.  The  first  one  is  the  essentially  non-oscillatory  (ENO)  polynomial  inter¬ 
polation,  and  the  second  is  the  filtering  technique  for  Fourier  approximations. 

ENO  Polynomial  Interpolation 

We  will  follow  the  notation  used  in  [9].  Let  u(x)  be  a  function  defined  on  I  =  [0, 27r] 
and  {xi}^0  be  the  uniform  mesh  on  I,  X{  =  ih,h  =  j*-.  For  simplicity  of  illustration,  we 
assume  that  u(x)  has  only  one  discontinuity  at  x,  6  /.  Now  given  u(xt),  0  <  i  <  N,  define  a 
piecewise  m-th  order  polynomial  interpolant  Qm(x;  u)  for  u(x)  at  mesh  points  x;,  0  <  i  <  N 
as  follows: 

Qm(xi-,u)  —  u(x,)  for  0  <  i  <  N,  (2.6) 

and 

Qm(x iu)  =  qmJ+i(x,u)  for  Xj  <  x  <  Xj+i,  (2-7) 
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where  gmj+i(x;u)  is  a  polynomial  of  degree  m  defined  below. 

Polynomial  ?mj+i(x;u)  interpolates  u(x)  at  (m  +  1)  successive  points  x,,  im(j )  <  i  < 
im(j)+n t-  The  stencil  of  these  (m+1)  mesh  points  will  be  chosen  according  to  the  smoothness 
of  the  data  u(xi)  around  Xj.  A  recursive  algorithm  to  define  im(j )  starts  by  defining 


ii(j)  =  3  (2.8) 

i.e.  5lj+i(x)  will  be  the  first  degree  polynomial  which  interpolates  u(x)  at  Xj,xJ+1.  If  we 
assume  qkj+i(x)  is  the  k  -  th  degree  polynomial  which  interpolates  u(x)  at 

x*k(j)> ' '  ‘ )  x*k(j)+k>  (2-9) 


then  we  need  one  additional  mesh  point  in  order  to  define  qk+1j+\,(x).  That  point  may  be 
the  nearest  one  to  the  left  of  stencil  of  (2.9)  (i.e.  Xtk(j)-i)  or  the  nearest  one  to  the  right  of 
the  stencil  of  (2.9)  (i.e.  ^(jj+fc+i)-  The  choice  will  be  based  on  the  absolute  values  of  the 
corresponding  (A:  +  l)-th  order  divided  differences,  namely 


ik+i(j)  = 


HO')  1  if  U[x*fc(i)-1>  '  '  '  >  X‘)i0)+<:]  <  U[S«j.(j)>  ‘  '  '  i  x»j,(j')+A:  +  l]> 

ffc(i)  otherwise. 


(2.10) 


The  piecewise  polynomial  Qm(x]  u )  defined  in  (2.6),  (2.7)  will  give  uniform  nonoscillatory 
approximations  to  u(x)  up  to  the  discontinuities.  In  fact  it  can  be  shown  that 


-r-rQm{x]u)  =  —ru(x)  +  0(hTn+1  *)  for  0  <  k  <  m.  (2-11) 

dxk  dxk 

except  for  the  cell  containing  x,.  Using  sub-cell  resolutions  [10]  this  is  also  correct  in  the 
shocked  cell. 


Filtering  Techniques  for  Fourier  Approximations 

When  a  function  ti(x)  is  discontinuous,  its  Fourier  approximation  %(i)  will  be  at  most 
first  order  everywhere  [7]  .  However,  there  are  several  ways  to  recover  the  loss  of  the  spectral 
accuracy  in  the  smooth  regions  of  the  function  u(x)  ([8],  [14]).  The  most  common  way  is 
to  multiply  the  Fourier  coefficients  of  u(x)  by  a  decreasing  scalar  factor  a*.  |a*|  — >  0  as 
|fc|  — ►  N.  Then  the  resulting  series  (the  filtered  approximation)  will  be  denoted  by  u^(x), 


u£(x)  =  £  akake'kx.  (2.12) 

\k\<N 

It  was  proven  in  [20]  that,  if  ak  is  derived  from  a  scalar  function  cr(t),0  <  t  <  1  and 
=  °’('w  )>  1^1  —  -^i  an(f  CT(0  satisfies  the  followiiig  conditions: 
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(2.13) 


a(0)  =  1, 

a(l)  =  0, 

a<*\0)  =  a<fc>(l)  =  0  for  1  <  k  <  K, 

then  Utf(x)  will  converge  to  u(x)  in  the  smooth  regions  of  u(x)  in  the  order  of  0( jnrrr)- 
In  the  actual  computations,  ak  is  chosen  to  decay  exponentially  in  terms  of  the  frequency 
number, 


ak  =  for  |Jb|  <  N,  (2.14) 

where  the  constant  a  is  chosen  so  that  <7#  is  the  machine  zero  and  21  is  called  the  order  of 
the  exponential  filtering. 

Uniform  Spectral  Approximations 

Now  we  present  our  new  uniform  high  order  nonosciilatory  spectral  approximations  to 
discontinuous  functions.  Again,  for  simplicity,  we  assume  u(x)  is  a  periodic  piecewise  C°° 
function  on  [0,27r]  with  only  one  discontinuity  at  x,.  Also  we  assume  that  the  discontinuity 
has  been  detected  within  an  interval  [x*,x']. 

Let  us  denote  all  the  mesh  points  inside  interval  [x1,,  xr$\  as  Xj,,  •  •  • ,  XiT.  Then  we  define  a 
piecewise  m  -  th  order  polynomial  <p(x)  which  interpolates  the  function  u(x)  at  mesh  points 


'  qmj+  $(*)  if  x  e  lx}>  *i+i] n  [*i»  *;] for  some  h 

<p(x)  =  <  p,(x)  if  X  6  [0,x*],  (2.15) 

.  PT(x)  if  x  6  [xT„2i r], 

where  j+i(x)  was  defined  in  (2.8)  -  (2.10)  and  Pj(x )  and  PT{x)  are  both  m'-th  order 
polynomials  on  the  interval  [0,x*]  and  [x^,27r]  respectively  ,  m!  =  2m  +  1,  and  satisfy  the 
following  conditions, 


p[k\x\  -  2*) 

Pik\x\  +  2*) 


(2.17) 
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Conditions  (2.16)  (2.17)  ensure  that  yj(x)  will  be  at  least  globally  Cm  continuous.  There 
are  exact  2m  +  2  =  m'  +  1  constraints  on  the  m'-th  polynomials  Pi(x)  and  Pr(x )  respectively. 
Therefore  they  are  uniquely  defined.  By  (2.U),  the  function  <p(x)  will  have  the  following 
property, 


<p(x<)  =  u(xi)  for  ii  <  i  <  ir,  (2-18) 

<p(x)  —  u(x)  =  0(hm+1 )  for  x  E  [x*,x'].  (2.19) 

Next  we  consider  the  difference  between  u(x)  and  <p(x),  v(x)  =  u(x)  —  <p(x).  v(x)  will 
be  a  Cm  function  everywhere  in  [0,27t]  except  at  x,  where  [v(x)]Xj  =  0(/im+1).  Moreover, 
v(xi)  =  0,  ii  <  i  <  iT.  Therefore,  the  filtered  Fourier  interpolant  I^v(x)  will  converge  to 
v(x)  rapidly, 

T-1 

vn(x)  =  Inv(x)  =  akVkexkx, 

k=~r 

h  =  4l> 
iV  »=0 

=  4  EM*0  -  «(*())«•““  +  4  EW*i)  - 

and  <Tfc  is  the  filter  in  (2.14). 

Finally  we  define  the  uniform  spectral  approximation  Vu(x)  of  u(x)  by 

Vu(x)  —  <p(x)  +  v^(x),  for  x  G  [0, 27r].  (2.20) 

Then  the  derivatives  of  u(x)  will  be  approximated  by  those  of  Vu(x ),  i.e. 

~  Vu f°r  >  °'  (2'21) 


To  see  the  accuracy  of  (2.20)  to  u(x),  let  p(x)  E  C™(xlt,xTl)  be  a  mollifier  function  such 
that 


II 

TT 

for  x  near  x,, 

(2.22) 

and 

if  x  €  [xi.xj], 

otherwise. 

„.(I)  =  { :w»w 

(2.23) 

One  can  easily  see  that 

u*(x,)  =  v(xt) 

for  0  <  i  <  N, 

(2.24) 
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and  hence 


/jvu*(x)  =  Inv(  x).  (2.25) 

As  v*(x)  €  Cm(0,27r)  and  is  periodic,  by  standard  estimate  it  can  be  shown  [2]  that 


v-(x)  -  iwz)  Hi,  < 


iVm+1  ’ 


where  c  is  a  constant  independent  of  N . 
On  the  other  hand, 


by  (2.19)  we  have 


|u*(x)  -  v(x)||Ll  =  y ~  p(x))2v2(x)dx 

<  ^Jx\\x)dx 

<  Jjx‘(u(x)  -  v(x)y  dx, 


|w*(x)-w(x)||lj  =  0(/im+l). 


(2.26) 


(2.27) 


(2.28) 


It  follows  from  (2.26)  and  (2.28)  that 

IK*)  “  7W*)IU3  =  IK*)  -  ^v*(x)||ij 

<  IK*)  -  v'(x)\\L,  +  ||v*(x)  -  IfjV*(x)\\Li  (2.29) 

=  0(hm+1). 


Thus 


\\Vu(x)  -  titl)!!!,  =  ||[ip(l)  +  /Jt)(x)]  -  [^(l)  +  v(l)]j|tl 

=  IK*)  -  flU*)lk  =  0(/."*+'),  (2,30) 

which  establishes  the  uniform  (m+  l)-th  order  convergence  of  the  approximation  Vu  of  u  in 
the  Li  norm.  Error  estimates  in  a  higher  order  Sobolev  norm  can  be  derived  similarly. 

The  spectral  convergence  of  Vu(x)  to  tt(x)  in  the  regions  outside  [ x\,  x\ ]  follows  from  the 
spectral  convergence  of  v^{x)  to  v(x)  in  the  smooth  regions  of  v(x). 

3  Uniform  High  Order  Spectral  Methods 

In  this  section  we  study  uniform  high  order  spectral  methods  for  conservations  laws  (1.1). 
First,  we  will  consider  the  scalar  one-dimensional  conservation  laws.  Extensions  to  the 
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system  of  conservation  laws  and  to  multi-dimensional  problems  will  be  discussed  in  Section 

4. 

We  will  derive  the  spectral  schemes  using  the  method  of  lines.  The  time  derivative  and 
spatial  derivatives  will  be  discretized  separately.  For  simplicity  we  only  present  the  Euler- 
forward  difference  method  for  the  time  derivative.  In  the  numerical  experiments  of  Section  5 
a  high  order  TVD  type  Runge-Kutta  time  discretization  [17]  is  used.  The  numerical  scheme 
will  be  written  in  the  following  conservative  form, 


,.n+l  _  n 


=  - -M/wj  - /i- i). 


(3.1) 


where  u"  ~  u(xj,  tn),Xj  =  jAx,tn  =  nAt,  and  Ax  and  At  are  the  spatial  mesh  size  and  the 
time  step  respectively,  A  =  /J+ 1  are  the  numerical  fluxes  . 

It  is  observed  [17]  that  if  there  is  a  function  h(x)  such  that 


then 


r,  ,  „  '>(•',+  !)  -  *(**-*) 

fMx,))  = - ^ - ■ 


(3.2) 


(3.3) 


This  suggests  that  the  numerical  flux  f-+ 1  should  approximate  h(x]+^)  as  Ax  — >  0. 

We  construct  h(x)  in  the  same  manner  as  in  [17]  via  its  primitive  function  H(x )  modulo 
a  linear  function, 

H{x)  =  (3.4) 


where  c  is  a  constant  chosen  so  that  H(x)  will  be  a  periodic  function: 


-i: 


Ae 

T 


h(()d(  =  Ax^fr 
j= o 


(3.5) 


Assuming  that  (3.2)  holds,  then 

«(Vj)  =  /'”*(*(  ()-')<% 

J  3 


(3.6) 


] 

=  Ax  ^  /(u(xfc))  —  c(j  4-  l)Ax  for  0  <  j  <  N. 
k  =  o 

We  then  form  the  uniform  spectral  approximation  operator  VH  to  H(x), 

VH  =  v?(x)  +  u*-(x), 


(3.7) 
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where  <^>(x)  is  the  piecewise  m-th  polynomial  defined  in  (2.15)  and  u)y(x)  is  the  filtered 
Fourier  interpolant  of  data  H(xJ+i )  —  <p(xj+±)i  0  <  j  <  N-  As  before,  it  is  assumed  that 
the  shock  discontinuity  or  contact  discontinuity  x,  (for  simplicity  of  illustration,  only  one 
such  discontinuity  is  assumed  to  exist)  has  been  detected  in  an  interval  [x[,xra],  i.e. 

x»  €  (3-8) 

If  x  G  [Xj,xj+1]  n  [x',z;]  for  some  j,  <p(x)  =  qmj+  i(x).  gmii+i( x)  is  a  m-th  order 
polynomial  and 


i for  iJl)  <  i  <  i„(j)  +  m. 


(3.9) 


The  stencil  ■  ■  ■  ,im{j )  +  m  is  defined  recursively  as  in  (2.8)  -  (2.10),  however  the  first 

point  of  the  stencil  i\(j)  is  chosen  according  to  the  local  Roe-speed  a;+ i , 


i.e. 


*i+ 


1 

2 


f{uj- fi)  ~  fM 

Ui+1  -  Uj 


(3.10) 


if  ai+i 


if  a. 


J+5 


>0, 

<  0. 


(3.11) 


Then  we  have  the  following  spectral  algorithm 


Algorithm  I  (Spectral  ENO-Roe) 

t  Step  1,  define  f7(xJ+ 1),0  <  j  <  N  by  (3.6)  and  their  uniform  spectral  approximation 
VH(x)  by  (3.7); 

•  Step  2,  let 

fj+l  =  +  C  =  +  ^V^xi+\)  +  c>  (3’12) 

where  constant  c  is  defined  in  (3.5). 

Remarks 

1.  The  evaluation  of  £<«ii+i),0  <]  <  N  can  be  performed  efficiently  via  the 
Fast  Fourier  Transforms  with  the  total  number  of  operations  of  order  0(N  log  AT); 

2.  Regarding  the  first  term,  ^<£>(x^+i ),  o<;  <  N,  only  for  those  xJ+i  G  [x\,  xTa],  do 
we  need  to  do  the  ENO  interpolation  and  the  numerical  differentiation.  If  xJ+i 
is  located  outside  [x!f,x^],  the  derivatives  can  be  computed  analytically  without 
additional  computer  cost; 
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3.  The  formal  spatial  accuracy  of  Algorithm  I  will  be  spectrally  accurate  in  the 
smooth  regions  of  the  solution  and  uniformly  m-order  elsewhere; 

4.  (Entropy  Correction)  the  fluxes  /J+i  defined  in  (3.12)  are  based  on  the  Roe  flux 
which  admits  “expansion  shocks”.  We  use  an  entropy  correction  discussed  in 
[17],  see  also  Harten  [11]  for  the  origination  of  such  corrections. 

4  Euler  Equations  of  Gas  Dynamics 

In  this  section  we  extend  the  scalar  Algorithm  I  from  the  previous  section  to  the  system  of 
Euler  equations  for  gas  dynamics  for  polytropic  gas.  With  all  variables  in  bold  face  denoting 
vectors,  we  have  the  following  Euler  equations: 


u  t  +  f(u)*  =  0,  (4.1) 

u  =  (p,  m,  E)T ,  (4.2) 

f(u)  =  qu  +  (0,P,qP)T,  (4.3) 

P  =  (7  -  1)(-S  -  |p92)-  (4-4) 


Here  p,  q,  P  and  E  are  the  density,  velocity,  pressure  and  total  energy  respectively,  m  =  pq 
is  the  momentum  and  7  =  1.4  is  the  ratio  of  specific  heats  of  a  polytropic  gas. 

The  eigenvalues  of  the  Jacobian  matrix  A(u)  =  f£  are 


Mu)  =  q~c,  Aa(u)  =  q,  A3(u)  =  q  +  c, 


(4.5) 


where  c  =  (7 P/p)*  is  the  sound  speed. 

The  corresponding  right-eigenvectors  are 


rl(u)  = 


(  1  ^ 

(  1  ^ 

f  1  \ 

q-c 

,r2(u)  = 

,r3(u)  = 

q  +  c 

V  H  —  qc  ) 

U?2 ) 

\  H  +  qc  ) 

(4.6) 


where 


H  =  (E  +  P)/p  =  €*/(■  7-1)  +  !?’, 

is  the  enthalpy. 

The  corresponding  left  eigenvectors  {lfc(u)}  which  are  bi- orthonormal  to  {r*(u)}  in  (4.6) 
are 


Mu)  =  +  g/c,-6i9  -  1/c,  61), 

Mu)  -  (1-MM,-M>  (4-7) 

Mu)  =  ^(*2  -9/c,-M  +  1/c,&i)> 
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where 


(4.8) 

(4.9) 


1>1  =  (tt  -  1  )/c1, 

b2  1 

As  in  the  case  of  scalar  conservation  laws,  we  first  define  a  vector  counterpart  H(xy+i), 
0  <  j <  N  of  (3.6).  The  scalar  quantities  /(•u(xfc))  in  (3.6)  are  replaced  by  the  vectors 
f(u(xjt)).  We  then  generalize  the  construction  of  the  uniform  spectral  approximation  oper¬ 
ator  VH  to  vector- valued  functions  in  the  following  way.  As  before,  the  assumption  about 
the  shock  location  (3.8)  still  holds  here.  We  have  the  uniform  spectral  approximation 

?H(x)  =  $(x)  +  v^),  (4.10) 

where  the  components  of  the  vector-valued  function  $(x)  will  be  piecewise  m-th  order  poly¬ 
nomials  and  vjy(x)  will  be  the  filtered  Fourier  interpolant  of  vector  quantities  H(xJ+i)  — 

$(2j+i),0  <  j  <  N. 

3>(x)  will  be  defined  separately  according  to  whether  x  is  inside  or  outside  the  interval 
[xla,xr,].  For  x  £  [xj,xJ+i]  D  [x',x'j  for  some  j,  first  let  q^-+i(x)  be  the  ENO  polynomial 
interpolant  of  the  characteristic  variables  H.(k\xi+i),  j  —  m  <  i  <  j  +  m.  Here  as  usual 
the  characteristic  variables  H^fe^(xi+i)  are  the  projections  of  H(xi+^)  on  the  locally  defined 
characteristic  fields.  In  this  paper,  we  define  the  local  characteristic  fields  with  respect  to 
the  Roe-averaged  state  uJ+i  between  the  states  Uy,  Uy+i,  i.e. 

Hl‘>(x1+i)  =  U(u,+  ,)H(xj+i)  (4.11) 

for  j  —  m  <  i  <  j  +  m,  1  <  k  <  3.  See  [15]  for  the  definition  of  Uy+i. 

We  thus  have 


=  H(fc)(x»+* )  for  im{j)  <  i  <  im(j)  +  m.  (4.12) 

and,  in  the  recursive  process  of  choosing  the  stencil  im(j),  the  first  point  ii(j)  is  determined 
according  to  the  sign  of  the  local  eigenvalue: 


if  AW(Uy+l)  >  0, 
if  AW(u;+i)  <  0. 


(4.13) 


We  then  define  the  vector- valued  function  $(x)  as 


*(*)  =  ^j+i(x)n(uJ+i)  for  x  G  [xy,  Xj+J]. 
*=i  ’  3 


(4.14) 
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On  the  other  hand,  when  x  is  outside  of  the  interval  [x*,x'],  we  define  $(x)  in  the  same 
way  as  in  (2.16)  and  (2.17).  Therefore,  in  those  regions  the  components  of  $(x)  will  be 
m'  =  2m  -f  1-th  order  polynomials.  Globally,  the  components  of  $(x)  will  be  Cm  function 
for  x  €  [0,  2tt]  and 

*(*<+*)  =  H(*i+ i).  if  *  G  [*'.,*;] •  (4.15) 

Now  we  present  our  algorithm  for  (4.1) 
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Algorithm  II  (Spectral-ENO-Roe) 


•  Step  1,  define  vector  quantities  H(xy+i),  0  <  j  <  N  by  (3.6)  and  their  uniform  spectral 
approximation  VH(x)  by  (4.10); 


•  Step  2,  let 

U  =  *)  +  c  +  s^»)  +  * 

where  c  is  defined  as  in  (3.5)  with  fj  replaced  by  fy. 


(4.16) 


Remarks 

1.  All  of  the  remarks  following  Algorithm  I  of  the  previous  section  apply  to  Algo¬ 
rithm  II.  However,  the  entropy  correction  is  used  only  on  the  genuine  nonlinear 
fields.  The  total  number  of  operations  in  the  computations  of  all  f -+i ,  0  <  j  <  N 
will  be  of  order  0(N  log  N).  Moreover,  the  characteristic  decompositions  are 
only  needed  for  those  xJ+i  in  [xls,xTe]  where  a  shock  discontinuity  or  contact 
discontinuity  has  been  detected; 

2.  (Generalization  to  Mulit- dimensional  Cases)  For  the  system  of  conservation  laws 
in  two  dimensions  (1.1),  we  apply  Algorithm  II  to  f(u)  and  g(u)  separately. 
Characteristic  decompositions  will  be  performed  on  their  corresponding  charac¬ 
teristic  fields  when  needed.  The  same  idea  can  be  applied  to  higher  dimensional 
cases. 


5  Numerical  Results 

In  this  section,  we  will  carry  out  several  numerical  experiments  with  Algorithm  I  &  II.  In 
implementing  these  algorithms,  we  have  to  choose  the  order  of  ENO  interpolations  in  (3.7) 
or  (4.10)  and  the  interval  [ xls,xra j  as  defined  in  (3.8),  which  is  detected  to  contain  a  shock 
or  contact  discontinuity.  In  most  of  the  tests,  we  choose  third  order  ENO  interpolation,  i.e. 
m  —  3,  unless  it  is  mentioned  otherwise.  The  interval  [x1,,  x^]  is  usually  chosen  to  contain  6 
-  8  mesh  points  around  a  discontinuity.  The  numerical  results  show  insensitivity  to  the  size 
of  [x*,x^]  as  long  as  it  contains  all  the  transition  points  in  the  numerical  shock.  To  detect 
the  shock  we  have  used  the  basic  check  on  the  gradients  of  the  numerical  data.  Define 

tj  =  max(\u(j)  -  u(j  -  1)|,  |u(;  +  1)  -  u(j)|)  (5.1) 
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If  tj  >  max(3.0ij_2, 3.0tj+2,  a),  where  a  is  chosen  dynamically  according  to  the  structure 
of  the  shock  wave,  then  we  decide  that  the  interval  [xj-i,Xj+i]  contains  a  discontinuity. 
Because  of  the  global  high  order  of  the  scheme,  a  falsely  detected  shock  location  in  a  smooth 
region  of  the  solution  would  not  destroy  the  overall  accuracy  of  the  scheme. 

To  retain  the  spectral  accuracy  in  the  smooth  region  of  the  solutions,  we  apply  high  order 
exponential  cut-off  filters  in  (2.14).  It  is  our  experience  that  a  very  weak  filter  (i.e.  high 
order)  will  suffice  to  get  high  accuracy  in  the  smooth  region. 

Time  Discretization  for  Chebyshev  Methods 


The  time  derivative  in  (4.1)  is  discretized  with  the  Runge-Kutta  method.  We  have  used 
the  third  order  Runge-Kutta  method  proposed  in  [17]  which  yields  TVD  (total  variation 
vanishing)  results  if  the  spatial  discretization  is  TVD. 

For  periodical  problems,  the  spatial  derivative  is  approximated  by  Fourier  trigonometric 
polynomials.  When  the  solution  is  nonperodic,  Chebyshev  polynomials  are  used  instead.  A 
common  difficulty,  however,  with  Chebyshev  methods  is  the  stringent  time  step  restriction. 
In  general,  the  time  step  At  has  to  be  in  the  order  of  O(jp)  where  N  is  the  order  of  the 
Chebyshev  polynomials.  This  is  a  direct  consequence  of  the  clustering  of  Chebyshev  collo¬ 
cation  points  near  both  boundaries.  For  many  hyperbolic  problems,  this  dense  distribution 
of  mesh  points  near  boundaries  is  not  necessary,  as  for  most  of  the  test  problems  in  this 
paper.  Recently  in  [12],  a  novel  mesh  transformation  is  proposed  to  relax  the  restriction 
of  the  Chebyshev  methods  on  time  steps.  If  x  denotes  the  physical  coordinate  and  £  the 
computational  coordinate,  the  following  transformation  is  considered  in  [12]  : 


t  sin  1  ax 

£  ~  :  Tx 

sin  a 


If  &  =  cosijj,0  <  i  <  N  is  the  Chebyshev  mesh  in  the  space,  then  x,  =  isin(sin_1  a &) 
will  be  the  corresponding  mesh  in  the  x  — space.  Because  of  the  stretching  nature  of  the 
transformation  (5.2),  mesh  points  x*  will  be  more  uniformly  distributed  in  the  physical  x- 
space.  A  Chebyshev  polynomial  in  the  transformed  (-  space  will  be  used  to  approximate  the 
derivative  with  respect  to  £,  and  the  derivative  with  respect  to  x  will  be  computed  as  follows 


d  d  d£  a  d 

dx  d(  dx  sin-1  acos(sin_1  a£)  d£ 


(5.3) 


In  our  computations,  we  have  observed  an  improvement  of  one  magnitude  in  the  time  step 
with  a  =  0.999;  at  the  same  time,  the  resolution  of  the  numerical  method  is  also  enhanced 
in  the  interior  of  the  physical  domain.  We  refer  the  reader  to  [12]  for  more  details  about  the 
evaluation  of  the  transformation. 
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Uniform  Spectral  Approximation  for  Discontinuous  functions 

We  consider  the  following  piecewise  C°°  function 

—  sin(2(x  -f  0.77t))  +  1  —  tt  <  x  <  —  |7r, 

g(x)  =  e,mJ*  sin2  x  |x|  <  §tt,  (5.4) 

-  1  -  “  <  X  <  TT. 

^(s)  is  defined  periodically  on  the  outside  of  [ — 7r,  7r] ,  thus  g(x )  has  three  discontinuities  in 
the  interval  [0,27t].  The  profile  of  g(x )  is  plotted  as  the  solid  line  in  Fig. 2.  Now,  given  the 
mesh  value  g(xi),  X{  =  0  <  i  <  N,  we  use  (2.20)  and  (2.21)  to  approximate  g(x]+ 1) 

and  0  <  j  <  N  respectively.  In  Fig.  1(a),  (b),  we  plot  the  errors  in  function 

values  and  first  derivative  values  respectively  in  the  logarithm  scale  for  different  N’s,  i.e. 
N  =  32,64, 128,256.  In  all  of  the  runs,  we  use  a  16th  order  of  exponential  cut-off  filter  and 
the  third  order  of  ENO  interpolation  in  (2.15).  We  clearly  see  the  uniform  convergence  and 
spectral  accuracy  in  the  smooth  parts  of  the  function  g(x). 

Linear  Advection  of  Discontinuous  Solution  with  Subcell  Resolution 

In  order  to  reduce  the  smearing  in  contact  discontinuities  by  shock  capturing  schemes, 
Harten  [10]  suggested  a  sub-cell  resolution  technique  to  treat  one-dimensional  contact  dis¬ 
continuities  in  the  context  of  the  cell-averaged  ENO  finite  difference  schemes.  Later  on,  this 
idea  was  extended  to  the  point-value  version  of  the  ENO  finite  difference  scheme  in  [17].  We 
test  the  subcell  resolution  by  our  spectral  methods. 

Consider  the  initial  boundary  value  problem  of  the  following  linear  hyperbolic  equation 

Uf  —  ux, 

1  u(x,0)  =  g{x),  x  G  [0,27r]  (5.5) 

u(0,  t )  =  u(2n,  t ) 

where  g(x)  is  defined  in  (5.4). 

In  Fig.  2(a),  (b),  we  plot  the  numerical  solution  for  N  =  200  after  one  and  two  period¬ 
icities,  i.e.  t  =  27 r,t  =  47r.  In  both  runs,  we  have  used  the  10th  order  exponential  cut-off 
filters. 

Inviscid  Burgers  ’  Equation 

In  this  classic  example  of  shock  wave  computation,  we  consider  the  initial  value  problem 
with  sine  wave  initial  conditions, 

'  >*. + d). = o, 

•  u(x,  0)  =  a  4-  (3  sin  x  x  G  [0,27r],  (5.6) 

.  «(0,0  =  u(2n,t), 
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where  a  —  0.3,  (3  =  0.7. 

The  solution  to  (5.6)  develops  a  shock  discontinuity  at  time  t  =  jj.  When  a  ^  0,  the 
solution  consists  of  a  moving  shock  wave  after  t  =  L.  As  the  exact  solution  for  this  problem 
can  be  obtained  by  iterative  methods,  we  use  this  example  to  test  the  global  accuracy  of 
the  scalar  Algorithm  I  with  m  =  3  (i.e.  the  third  order  of  ENO  interpolation  is  used  in 
(3.7)).  In  table  1,  we  list  the  global  L\  error  of  the  numerical  solutions  and  the  L\  error 
in  the  smooth  regior  of  the  solution  at  time  t  =  2  and  N  =  32,64,128,256.  In  computing 
the  global  L i  error,  we  exclude  the  occasional  one  transition  point  across  the  shock,  and 
the  smooth  region  is  taken  as  the  region  which  is  0.8  away  from  the  shock  location.  The 
third  column  of  Table  1  shows  the  global  third  order  accuracy  in  L\  norm  of  Algorithm 
I,  and  the  fifth  column  shows  the  increasing  order  of  accuracy  in  the  smooth  region.  In  all 
of  the  runs,  the  time  step  has  been  chosen  such  that  further  decreasing  of  the  time  step 
does  not  improve  the  final  accuracy.  Therefore,  the  dominant  errors  come  from  the  spatial 
discretization.  In  Fig.  3,  solutions  at  time  t  =  2  and  N  =  32  are  plotted  (plus)  against  the 
exact  solution(solid  lines).  In  Fig.  4(a),  the  errors  for  N  =  32,  64, 128, 256  at  time  t  =  2  are 
plotted  in  logarithm  scale.  For  a  comparison,  we  plot  the  errors  with  the  same  parameters 
for  the  third  order  ENO  finite  difference  scheme  in  Fig.  4(b). 


One- Dimensional  Euler  Equations  for  Gas  Dynamics 

We  solve  the  system  of  equations  (4.1)  with  different  initial  data. 

1)  Interaction  of  a  shock  wave  and  density  waves 

We  consider  the  following  initial  condition  for  (4.1), 

(  Pi  =  3.857143,  qi  =  2.629367,  Pi  =  10.333333  -1  <  x  <  -0.8,  ,  * 

|  pT  —  1  +  esin57rx,  qr  =  0,  PT  —  1  —0.8  <  x  <  1, 

where  e  =  0.2. 

In  Algorithm  II,  we  choose  the  third  order  ENO  interpolation  and  the  16th  order 
exponential  cut-off  filters.  In  Fig.  5(a),  we  display  the  density  profile  for  N  =  200  at 
t  =  0.36,  the  solid  lines  are  the  solutions  obtained  by  the  third  order  ENO  finite  difference 
methods  with  800  points.  For  a  comparison,  we  plot  the  solution  obtained  with  the  second 
order  MUSCL  schemes  with  N  =  200  in  Fig.  5(b). 


2)  Sod’s  problem  and  Lax’s  problem 

We  now  consider  the  standard  Riemann  problem  of  (4.1)  with  the  following  initial  data 
[16], 

a)  Sod’s  problem 
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b)  Lax’s  Problem 


(. PhQhPi )  =  (1,0,1)  -l<z<0, 

(/ pr,qT,PT )  =  (0.125,0,0.10)  0  <  x  <  1; 


(. pi,qi,Pi )  =  (0.445,0.698,3.528)  -  1  <  x  <  0, 

(, Pr,qr,Pr )  =  (0.5,0,0.571)  0  <  x  <  1. 

For  the  Lax’s  problem,  the  sub-cell  resolution  has  been  applied  on  the  linear  degenerated 
field  near  contact  discontinuities.  In  both  of  the  problems,  we  use  the  third  order  ENO 
interpolation  in  Algorithm  II  and  the  10th  order  exponential  cut-off  filters.  In  Fig.  6 
(a)-(c)  the  solutions  of  the  Sod’s  problems  with  N  =  150  at  time  t  =  0.4  are  plotted,  while 
Fig.  7(a)-(c)  display  the  solutions  of  the  Lax’s  problem  with  N  =  150  at  time  t  =  0.26.  In 
both  cases,  the  solid  lines  are  the  exact  solutions. 

Interaction  of  blast  waves 

The  initial  data  suggested  in  [20]  to  simulate  the  interactions  of  two  blast  waves  is  as 
follows 


u  L 

0  <  X  <  0.1, 

UAr 

0.1  <  x  <  0.9, 

(5,8) 

Ufl 

0.9  <  x  <  1, 

where  pL  =  pM  =  pR  =  1,  qL  =  qM  =  qR  -  0,  PL  =  103,  Pm  =  10  2,Pr  =  102. 

The  solutions  to  this  problem  possess  drastic  fluctuations  under  the  impact  of  inter¬ 
actions;  it  is  a  good  test  of  the  stability  of  Algorithm  II.  The  complex  structure  of  the 
solutions  after  the  clash  of  two  blast  waves  demands  a  stable  high  order  method  to  capture 
the  details  of  solutions.  Unlike  the  finite  difference  methods,  the  spectral  methods  do  not 
require  exterior  mesh  points  to  treat  boundary  conditions.  We  apply  characteristic  boundary 
conditions  on  both  boundaries.  As  the  boundaries  are  treated  as  solid  walls,  we  impose  the 
condition  that  velocity  variables  vanish  on  both  boundaries,  i.e.  qo  =  0 ,qN—  0. 

In  Fig.  8(a)-(c)  and  Fig.  9  (a)-(c),  we  plot  the  solutions  of  the  state  variables  with 
N  =  300  at  time  t  —  0.028  and  t  =  0.038,  respectively.  The  former  is  an  instant  before  the 
clash  and  the  latter  is  one  after  the  clash.  The  solid  lines  in  both  Fig.  8  and  Fig.  9  are  the 
solutions  obtained  with  the  third  order  ENO  finite  difference  methods  with  800  mesh  points. 
In  Fig.  10,  the  solutions  of  density  with  N  =  400  are  also  plotted. 

Interaction  between  a  Two-dimensional  Shock  Wave  and  a  Rotating  Vortex 
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The  equations  in  consideration  will  be  (1.1)  with 


u  =  ( p,mx,my,E ),  (5.9) 

f(u)  =  qxu  +  (0,  P,0,qxP),  (5.10) 

g(u)  =  9vu  +  (0, 0,  P,qyP),  (5.11) 

where  qx,qv  are  velocity  components  in  x-  and  y  -  directions  respectively,  mx  =  pqx  and 
mv  =  pqy  are  x  -  and  y  -  momentums  respectively  .  P  =  (7  —  l)(i£  —  \pq2)>  Q2  —  ql  +  9y- 


We  apply  one-dimensional  Algorithm  II  on  the  fluxes  f(u)  and  g(u)  separately.  The 
right  and  left  eigenvectors  for  the  Jacobian  matrices  Jj,  |*  can  be  found  in  [17]. 

The  physical  domain  will  be  the  rectangle  [0,3]  x  [—1.5, 1.5].  A  Mach  3  planar  shock 
wave  moves  from  the  left  to  the  right.  A  rotating  vortex  is  initially  located  to  the  right  side 
of  the  shock.  As  the  time  progresses,  the  shock  will  hit  the  vortex  and  interact  with  it.  The 
shock  front  will  be  deformed  by  the  interaction,  and  pressure  waves  are  generated  from  the 
interactions.  In  the  computations,  we  define  the  velocity  fields  of  the  vortex  as  those  induced 
by  two  rotating  concentric  cylinders  with  radius  ri  and  r2  respectively,  rx  <  r2.  Initially  the 
vortex  is  located  at  ( xc,yc ).  The  outside  cylinder  is  stationary  and  the  inside  one  rotates 
with  the  angular  velocity  1 0.  Let  t;(r)  be  the  radius  velocity  at  a  distance  r  from  the  center 
of  the  vortex,  we  have 

'  c or  if  0  <  r  <  ri, 

u(r)  =  «  WI(.L  IIl)  if  n  <  r  <  r2,  (5.12) 

.0  if  r  >  r2, 


where  a  =  b  =  r\  -  rf,  r  =  \J{x  -  xc )2  +  {y  -  yc )2.  We  choose  rj  =  0.15,  r2  = 

0.75, u;  =  7.5.1 

Therefore  the  x  -  and  y  -  velocities  induced  by  this  vortex  at  ( x,y )  will  be 


y  -  yc.,  s 
?*  =  ■y(T’), 

T 

(5.13) 

X  —  xc  _  f  ^ 

9v  =  v\r)i 

r 

(5.14) 

where  xc  =  2.25,  yc  -  0. 

The  initial  conditions  for  the  simulation  will  be  as  follows: 

pi  =  3.857143, 

qxi  =  2.629367  if  x  <  x0, 

€ 

11 

0 

(5.15) 

Pi  =  10.333333, 

(5.16) 
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and 


Pr  —  1  > 

9ir  * —  Qx  if  X  x*  Xo, 

*}yr  =  4y>  (5.17) 

Pr  =  1, 

where  x0  is  the  initial  shock  position,  Xo  =  1.5. 

We  impose  characteristic  boundary  conditions  on  both  left  and  right  boundaries.  A 
periodic  boundary  condition  is  considered  in  the  y-direction  and  hence  we  are  actually  sim¬ 
ulating  the  interaction  between  an  array  of  periodically  distributed  vortices  and  a  planner 
shock  wave.  To  relax  the  time  restrictions  of  the  Chebyshev  approximation  in  the  x  -  di¬ 
rection,  we  apply  the  mesh  transformation  (5.2)  with  a  =  0.999.  The  shock  has  been  made 
stationary  by  a  translation  in  the  mean  flow  direction.  The  second  order  ENO  interpolation 
and  the  10-th  order  exponential  filter  have  been  used  in  (4.10). 

Fig.  ll(a),(b)  are  the  contour  plots  of  the  pressure  and  density  fields  at  time  t  —  0.4, 
while  Fig.  12(a), (b)  are  the  close-ups  of  the  pressure  and  density  at  time  t  =  0.4.  Fig.  13  is 
the  pressure  profile  at  time  t  —  0.4. 

Concluding  Remarks 

Centered  difference  methods  including  spectral  methods  are  efficient  and  accurate,  while 
upwinding  difference  methods  offer  the  advantages  of  sharp  monotonic  shock  profiles.  We 
have  explored,  in  this  work,  the  possibilities  of  blending  the  advantages  of  the  ENO  finite 
difference  methods  and  the  spectral  methods.  Numerical  results  have  shown  the  robustness 
and  feasibility  of  this  approach  at  a  small  extra  cost  over  the  standard  spectral  methods. 
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Table  1:  Global  L\  error  and  Zq  error  in  the  smooth  region  for  the  Burgers’  equations  at 
time  t  =  2,  the  smooth  region  is  defined  to  be  0.8  away  from  the  shock  location 


N 

Global 

Order 

Smooth  region 

Order 

32 

64 

1.49(-4) 

2.5 

1.17(-4) 

4.3 

128 

2.70(-5) 

2.9 

5.86(-6) 

6.5 

256 

3.70(-6) 

3.7 

6.54(-8) 

10.0 

2.95(-7) 

6.36(- 11) 
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on  of  Discontinuous  Solutions  with  Subcell  Resolutions 


(a)  Density,  (b)  Velocic 


the  Lax's  Problem  with  N  =  150  at  time  t  =  0.26.  (a)  density,  (b)  velocity 


Figure  8.  (Continued) 
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Figure  8.  (Continued) 


A.  x  d_  s 


Figure  9.  Blast  Wave  Problem  with  N  =  300  at  Time  t  =  0.038.  (a)  Density,  (b)Velocity,  (c)  Pr 


s  Betwwen  a  Two-Dimensional  Planner  Shock  Wave  and  A  Vorcex,  Shock  Mach  numer 
y  Contour  with  Level  Value  =0.1,  (b)  Pressure  Contour  with  Level  Value  =  0.2 


Figure  12.  Close-ups 


48 


Figure  12.  (Conti 


as  in  Figure  LI,  Densicv  Solutions  at  Time 
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